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I. Introduction 

This paper investigates connections between fractional 
viscoelasticity, fractional wave equations, causal models, 
and power-law attenuation within the framework of elas- 
tic wave modeling. Paying extra attention to medical 
imaging applications, we intend to convey to the frac- 
tional calculus community information on developments 
related to time-fractional elastic wave equations. The 
paper expands upon the conference proceeding paper [l| 
and the J. Acoust. Soc. Am. papers 0, 

In Section ini the frequently encountered power-law at- 
tenuation nature of elastic waves in complex media is 
considered, especially in medical imaging. 

In Section lllli a fractional generalization of the Zener 
viscoelastic model is reviewed. Parameter value restric- 
tions to keep the model thermodynamically admissi- 
ble are discussed and experimental evidences are listed. 
The section also demonstrates the important achieve- 
ment that the fractional Zener model is equivalent to 
a Maxwell- Wiechert representation consisting of a set of 
conventional springs and dashpots. Also physically un- 
derlying principles which lead to fractional viscoelasticity 
are reflected upon. 

Subsequently, Section IIVI derives a fractional Zener 
wave equation from conservation laws and the fractional 
Zener constitutive stress-strain model. Section |V] then 
analyzes properties of this wave equation. A connec- 



tion to the widely acknowledged Nachman-Smith-Waag 
wave equation for acoustic propagation with relaxation 
losses [J] is demonstrated, attenuation and finite phase 
speed power-law regimes are evaluated, and causality con- 
ditions are considered. 

The final section provides conclusions and discussions 
including parallels to similar models in other fields, such 
as electromagnetics. 

II. Power-law attenuation in complex media 

Elastic wave attenuation in complex media such as bi- 
ological tissue, polymers, rocks, and rubber often follows 
a frequency power law: 

ak{io)o,io\ (2.1) 

with the exponent 77 e [0,2]. Such power-laws can be 
valid over many frequency decades and examples are 
found all the way from infrasound to ultrasound [5|. See 
e.g. Fig. 1 in [6] which visualizes experimentally estab- 
lished power-law attenuation examples for both shear 
and compressional waves. As summarized in [7], com- 
pressional wave attenuation in biological tissue commonly 
manifests a power- law exponent ry e [1,2], while for shear 
waves in tissue one often finds 77 G [0, 1]. 

A. Power-law attenuation in medical imaging 
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The establishment of accurate wave-propagation mod- 
els in power-law lossy media is important for applications 
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where broadband waves are utihzed. Below we hst exam- 
ples of such applications. 

In photoacoustic imaging and tomography, laser pulses 
are delivered into tissue where an associated thermoelas- 
tic expansion induces a wideband ultrasound emission 
which is used for image reconstruction [sj-flTj. 

Magnetic resonance (MR) elastography is a another 
means to estimate soft tissue stiffness e.g. in the liver or 
the breast. This method utilizes the propagation of shear 
waves which are monitored using MR imaging [l^ - [l8l |. 
Doppler techniques may also be applied to monitor the 
frequency-dependent tissue shear wave properties [Toj . 

A related elastography method is acoustic radiation 
force imaging (ARFI) techniques, where tissue is de- 
formed by a short compressional pulse which creates a 
propagating shear wave. The displacement due to the 
shear wave as measured by ultrasound is then used to 
quantify the tissue's mechanical properties e.g. by esti- 
mating the shear wave propagation speed |20l422| . 



B. Modeling power-law attenuation 

For acoustic modeling, time-fractional derivative wave 
equations have been shown to imply power-law attenu- 
ation over wide frequency bands [2, As further de- 
scribed in Section IIV[ fractional wave equations can be 
obtained from linearized conservations of mass and mo- 
mentum in combination with time-fractional constitutive 
relations connecting stress and strain. Related linear 
wave-propagation simulations are reported e.g. in p3l - l26| . 
For waves at high amplitudes, non-linear effects need to 
be considered. Such models were developed in [l^, [1^ 
where the fractional Kelvin- Voigt constitutive equation 
was applied. See also the related recent paper [29|. 

From a numerical modeling point of view when consid- 
ering the low-frequency or small-attenuation regimes of 
a power-law attenuating medium, wave equations with a 
d'Alambertian part and time-fractional terms may conve- 
niently be converted into space- fractional Laplacian mod- 
els. This can be beneficial due to reduced time-signal 
storage needed in the propagation simulation steps [29|- 
[3I I . Special care needs to be taken to ensure causality of 
the resulting models. The authors are not aware of simi- 
lar valid conversions between fractional temporal deriva- 
tives and fractional Laplacians that are applicable for at- 
tenuation with 7y < 1 or within a high-frequency regime. 

The multiple relaxation mechanism framework of Q is 
widely considered as adequate for acoustic wave modeling 
in lossy complex media like those encountered in medical 
ultrasound. It relies on thermodynamics and first princi- 
ples of acoustics. The corresponding wave equation for 
N relaxation mechanisms is a causal partial differential 
equation with its highest time derivative order iV-|-2. We 
denote this the Nachman-Smith-Waag (NSW) model. In 
order to make the discrete NSW model attenuation ade- 
quately follow w'', either the valid frequency band must 
be narrow, or the number of assumed mechanisms N 



must be large thus inferring a partial differential equa- 
tion of high order. 

On the other hand, for a certain continuous distribu- 
tion of relaxation mechanisms, the NSW and fractional 
Zener descriptions are equal and power-law attenuation 
is attained, as further reviewed in Section IV Al 

Band-limited fits to power-law acoustic attenuation for 
relaxation models with N = 2 and 3 are exemplified by 
[H,!!!]. In the latter, one of the mechanisms is assumed 
to be of very high relaxation frequency, thus represent- 
ing a thermoviscous component. The determination of 
the two other relaxation frequencies and their compress- 
ibilities, as well as the compressibility contribution of the 
thermoviscous component are then decided by numerical 
minimization of the resulting difference to the desired 
power-law attenuation. For a large number of modeled 
relaxation mechanisms, such numerical optimization of 
the parameter fit turns very intricate. 

Another approach is used in (35l |. which is closely re- 
lated to [H, H^l . It demonstrates that hierarchical frac- 
tal ladder networks of springs and dashpots can lead 
to power-law acoustical attenuation in a low-frequency 
regime. This however requires a large number of degrees 
of freedom which makes parameter fits cumbersome. 



III. The fractional Zener constitutive relation 

The history of fractional derivatives in physics goes 
back to Abel's integral equation from 1826 [38], which 
turns out to correspond to the 1/2-ordcr derivative. 
Early viscoelasticity- related papers are [39|, |40|, see also 
historical overviews in [4l|, Inclusion of fractional 

derivatives in the viscoelastic stress-strain relationship is 
convenient for describing many materials where the re- 
sponse depends on the past history, see e.g. [4l|, |43| 
(for reflections on this from an acoustical point of view, 
see [isl). For a record on the most intuitive fractional 
generalizations of the conventional (no n- fractional) vis- 
coelastic models, see e.g. the survey [40|- In addition, the 
comprehensive reviews (47l . |48| , summarize research on 
fractional calculus in dynamic problems of solid mechan- 
ics. Illustrations of the fractional derivative viscoelastic 
models commonly include the spring-pot (or just pot) el- 
ement. 



A. Stress-strain relation 

A five-parameter fractional Zener model fractional gen- 
eralization may be expressed as 



a{t) 



df^ajt) 



^^d^e{t) 



(3.2) 



where t is the time, a{t) the stress, e{t) the strain, r^r and 
Te positive time constants, and Eq the modulus. Here the 
nomenclature for the time-fractional derivatives follows 
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[49|, however many authors display naming conventions 
where a and (3 are interchanged. To be physically ade- 
quate, one requires a = {3 as further investigated in Sec- 
tion HTlB] and Section fV Al The fractional Kelvin- Voigt 
constitutive relation may be regarded as a low-and inter- 
mediate frequency representation of the fractional Zener 
model [2*1 which corresponds to — > in p.2p . 

Other fractional stress-strain relations with either the 
same or more degrees of freedom may be used to describe 
material response, as stated in e.g. |48|. One example is 
the 5-parameter approach described in (soj . Such and 
other generalized models could equally well be applied in 
the wave equation derivations described in the following. 

B. Parameter value restrictions 

Based on arguments from [s!], [l^l , a monotonically de- 
creasing stress relaxation requires a — fi \\\ the stress- 
strain relation p.2p . Table U lists thermodynamical of 
constraints on the (|3.2p parameters. See also the related 
recent paper fsslj. 

TABLE I: Fractional Zener stress-strain model p.2|) parame- 
ter constraints, from f5l!|. 

Eo>Q, EoT^>0, TS>r!, > Q, a = P 



We note that even though the rheological model (|3.2p 
is not thermodynamically well-behaved for a ^ fi^ the 
corresponding underlying mechanical model can actually 
be admissible. However this is only if the instantaneous 
wave is allowed to propagate at infinite speed [sij. 

The fractional Zener model with a = /3, which corre- 
sponds to the empirical Cole-Cole rheological relaxation 
spectrum formulation [sl] , was considered already in |39i] 
alongside with experimental evidences. 

C. Experimental evidences 

Parameter fits of experimental measurements to the 
fractional Zener model for biological materials include 
for brain [56H59[, human root dentin |60l|. cranial bone 
[6ll |. liver |57l|. arteries [sl], breast [63, and hamstring 
muscle l64|. Non-biological materials are exemplified by 
metals 1391, dop ed corning glass [i^, rubber [51], and 
polymers [6J, [65l - [69l | . For an account of experimental fits 
to fractional calculus stress-strain models made up to the 
year of 2009, see Section 2 in [4i|. 

The framework of viscoelasticity and acoustics in com- 
plex media has significant similarities with the framework 
of dielectrics and electromagnetic propagation. Not only 
do the wave equations have similarities in structure, but 
also the same constitutive relations connected to frac- 
tional derivative stress-strain modeling have experimen- 
tally been observed in the dielectrical properties of com- 
plex media. This is relevant for e.g. ground-penetrating 



radar and medical diagnosis using ultrawideband or Ter- 
ahertz electromagnetic waves. From an electromagnetic 
modeling point of view, the complex dielectric permittiv- 
ity plays the same role as the compressibility (the com- 
plex compliance) does in acoustics and viscoelasticity. 



D. Maxwell-Wiechert representation 

Viscoelastic constitutive stress-strain models are gen- 
erally possible to convert into a Prony series of Maxwell 
elements, that is a Maxwell-Wiechert description with 
springs and dashpots in parallel as illustrated in Fig. [T] 




FIG. 1: Illustration of the Maxwell-Wiechert viscoelastic rep- 
resentation. Damper i is denoted by Ei and dashpot i by r/i. 

Already in [tOI , the fractional Zener model was actually 
interpreted as a relaxation distribution. Several other au- 
thors have published related results where the fractional 
derivative stress-strain models are expressed without us- 
ing the exotic spring-pot viscoelastic element. In [7l| . a 
very large number of weighted Maxwell elements (Debye 
contributions) evenly distributed on a linear frequency 
scale are shown to give the same stress response as a 
fractional order viscoelastic model. Ref. [T^l, presented 
examples where viscoelastic damping due to several si- 
multaneously decaying processes with closely spaced ex- 
ponential decay rates are shown to induce a constitutive 
behavior involving fractional order derivatives. Further- 
more, Machado and Galhano have shown that averaging 
over a large population of microelements, each having 
integer-order nature, gives global dynamics with both in- 
teger and fractional dynamics i73|. We also note that 
the rheological fractional spring-pot element was inter- 
preted in terms of weighted springs and dashpots in [73 | . 
The recent paper [75| considers anomalous relaxation in 
dielectrics and interestingly provides relaxation distribu- 
tion functions for non-Debye models flavors such as the 
Cole-Cole, Davidson-Cole, and Havriliak-Negami, out 
of which the Cole-Cole one is most tightly connected to 
the current work because it corresponds to a — /3 in the 
fractional Zener viscoelastic model. 

As further discussed in Section IV A[ a Maxwell- 
Wiechert description of the fractional Zener model was 
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actually implicitly verified also in [s'] where it was con- 
nected to the multiple relaxation framework of [3| via a 
distribution with fractal properties. 



E. Physical background 

In [xl], the author reflects on the surprisingly good fit 
of the fractional Zener model to measured data, suggest- 
ing that this hints at the existence of underlying govern- 
ing principles. The proposed model is for internal energy 
loss at the molecular level of a viscoelastic polymer. A 
probability density to describe the motion of elevated en- 
ergy states along the molecule or "kinks" was postulated. 
Such states behave similar to particles in potential wells 
which have to overcome barriers. By certain assumptions 
on probability density functions for both an energy tran- 
sition and a distribution of barrier heights, it is shown 
that the probability density for an energy transition is 
fractal. The fractality leads to a relaxation function de- 
scribed by a power-law or a Mittag-Leffler function. This 
leads to a four parameter fractional Zener model linking 
stress and strain. 

Another interpretation is due to Mainardi who justifies 
the four-parameter fractional Zener model from the ther- 
moelastic equations when the temperature change due 
to diffusion and adiabatic strain change is a fractional 
derivative. The temperature change is thus a hidden vari- 
able in the stress-strain relationship fld\ . 

A third point of view considers propagating waves. At- 
tenuation of both compressional and shear waves is con- 
sidered to be due to two mechanisms: Absorption which 
is the conversion of energy into e.g. heat, and scattering 
which is the reflection of energy in all directions. Despite 
the different physical explanations, both mechanisms of- 
ten seem to result in power law attenuation. 

In medical ultrasound below about 10 MHz, it is gener- 
ally considered that absorption is the dominating mecha- 
nism [t^I. It is likely that the above models are relevant 
for this regime. For elastography in tissue with propa- 
gating shear waves in the 10 and 1000 Hz range, recent 
experimental results indicate that attenuation due to mul- 
tiple scattering can dominate [t^ [t^ . A 1-D model for 
multiple scattering attenuation is the O'Doherty-Anstey 
model [s^l, however there are no well established 2-D or 
3-D model equivalents. Assuming a fractal distribution of 
reflection coefficients, the O'Doherty-Anstey model can 
explain power- law attenuation (8l| . 

Connections between power-law attenuation, multiple 
scattering, and fractal geometry need to be further ex- 
plored in real 2-D and 3-D media. 



IV. Deriving the fractional wave equation 

Below we demonstrate that causal wave equations can 
be constructed from the expressions for linearized con- 
servation of momentum and the linearized conservation 



of mass combined with a fractional constitutive relation 
between stress and strain. The approach outlined was ap- 
plied already in [s^] to derive a wave equation based on 
the non-fractional Zener model, which is often denoted 
as the standard linear solid. 



A. Conservation laws 

In the following, we consider an isotropic medium 
where the only non-negligible stiffness parameters are ei- 
ther the bulk modulus, cn, or the shear modulus, C44, 
[S^. Then the linearized conservation of mass corre- 
sponds to the strain being defined by 



e(t) — \/u{x, t), < > e{uj) — ~ik u(k, w). 



(4.3) 



where u is the displacement in the transverse (for com- 
pressional waves) or longitudinal (for shear waves) direc- 
tions, X is the 3-D spatial coordinate, and the symbol 
J- denotes transformation into the spatio-temporal fre- 
quency domain where lo is the angular frequency and k 
the wave number. This is valid under the assumption of 
infinitesimal strains and rotations which is adequate e.g. 
in acoustical medical imaging. 

The linearized conservation of momentum is expressed 

as 



-ik(T{uj) ~ po{iuj)'^u{k,uj), 
(4.4) 



where po is the steady-state mass density and a denotes 
the stress, which in this context corresponds to the neg- 
ative of the pressure. 



B. Generalized compressibility k{uj) 

The frequency-domain generalized compressibility is 
defined as the ratio between strain and stress: k{uj) = 
e{uj)/a{uj), therefore being related to the constitutive 
stress-strain relation. Combining this definition with the 
conservation laws (14.31) and (14.41) gives 



V^u{x,t) 



n{t) * u{x^ t) < > k^{ijj) = uj^ P{)k{ijj). 

(4.5) 



Under circumstances where the linearized conserva- 
tions of mass and momentum are valid, the wave equa- 
tion is thus completely determined by the generalized 
compressibility. 

From p. 21) . the frequency-domain fractional Zener 
compressibility is obtained through the ratio e{u)) / a{uj): 



Kl{uj) = Kq 



1 -I- {T.ltjf 

1 -I- (r^iw)" 



= Kq — ILOKq- 



+ (jw)" 



(4.6) 
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The generalized compressibility k{uj) as given above 
is (e.g. in viscoelasticity) called complex compliance 
J*{uj) = 1/G*(w), where G*{uj) is the complex modulus. 

C. The fractional Zener wave equation 

Insertion of the generahzed compressibihty (j4.6p into 
the dispersion relation (14. 5p . results in the time-domain 
fractional Zener wave equation 

V. Properties of the fractional wave equation 

Below we explore some properties of the fractional 
Zener wave equation, first by connecting it to a mul- 
tiple relaxation wave model, then by studying the at- 
tenuation and the phase velocity as a function of fre- 
quency where 3 characteristic power-law regions are iden- 
tified. The causality of the model is finally verified. More 
mathematics-oriented studies of fractional Zener-related 
wave equations can be found in e.g. [84 - :86.] . 

A. Link to the NSW multiple relaxation framework 

Under the assumption of a continuum of relaxation 
mechanisms, the NSW model Q is linked to fractional 
derivative modeling in Q, where it was shown that the 
wave equation corresponding to a certain distribution 
of relaxation contributions is identical to the fractional 
Zener wave equation. As resumed in the following, the as- 
sociated compressibility contributions were shown to be 
distributed following a function related to the Mittag- 
Leffler function. The current section is conceptually 
tightly connected to Section IIII Dl where the fractional 
Zener constitutive model is represented as a set of springs 
and dashpots. 

1. The NSW generalized compressibility 

The NSW model of multiple discrete relaxation pro- 
cesses results in the generalized compressibility (which is 
equivalent to the rheological complex compliance J*{u})) 

N 

k{uj) = kq - iuS^ " , (5.8) 

u—l 

where the mechanisms v = 1 . . . iV, have the relaxation 
times Ti , . . ^ Tjv and the compressibility contributions 
Ki,...,KN [3. 

Note that (|5.8p corresponds to a sum of N weighted 
conventional Zener contributions, where each is given by 
dm with a = /3 = 1. 



Following fs*!, a representation of (|5.8I) when consid- 
ering a continuum of relaxation mechanisms distributed 
in the frequency band Q. G [f2i,f22] with the compress- 
ibility contributions described by the distribution K,y{Q) 
becomes 

K-f^ (lo) = kq — iu! [ — - dfl. (5-9) 
Joi " -f luj 

Letting the limits of the integral go between fii = 
and = 00, and instead incorporating any possible 
relaxation distribution bandwidth limitation into Ki^(ri), 
the integral above is a Stieltjes transform. Applying the 
Laplace transform relation 

/:.-{^}(i) = e-, (5.10) 

by virtue of Fubini's theorem [s^l the generalized com- 
pressibility (j5.9|) then becomes 

/•oo poo 

Jo Jo 
= Ko - iujTt[Hit)Cn {K,(17)Ki)}(^)- (5-11) 

Provided that the conservations of mass (|4.3p and mo- 
mentum (14. 4p are valid, and provided that the NSW 
generalized compressibility kn('^) of (I5.11|) is equal to 
the fractional Zener generalized compressibility kz(w) of 
(j4.6p . the dispersion relations from (j4.5p are also equal. 
Because the dispersion relation is a spatio-temporal 
Fourier representation of the wave equation, kn(w) = 
Kz(w) thus implies that the NSW wave equation becomes 
equal to the fractional Zener wave equation (14. 7p . Direct 
comparison of kn(w) in (|5.1ip to kz{uj) in (|4.6p . tells that 
they are equal in case the following is true: 

(za.)"-i - (rf /r,")(zc.)"-("-^+i) 



2. The Cole-Cole equivalent a — l3 case 

First we choose to study the case a = (3 in a similar 
manner as in . Inverse Fourier transformation of both 
sides of (|5.12p . then gives 

Hit)Cn {^Amit) = '^0(1 - r^/r:)^-' / ^^^^l^l!^ \ (t) 

L Ta + (iw)" J 

= Koil ~ T^/T^)H(t)E^^^ {-{tlr.r) , (5.13) 

where Ea^b{') is the Mittag-Leffler function (see Appendix 
IVII Al) . and Hit) is the Heaviside step function. The 
Fourier transform relation used in the last step above 
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is given in (j7.25p . Moreover, using the inverse Laplace 
transform relation of (j7.26L Eq. (j5.13p hence gives 



^o(l-r,7T,")/„,i(f7,T-") 
1 ko(t^ - sin(a7r) 
TT {rMy" + 2{tM)" cos(a7r) + 1 



(5.14) 



where fa,i{i^, a) was inserted from (j7.27l) . Note that this 
link between the fractional Zener and the NSW models is 
valid also outside the small-attenuation regime 9 {k} ^ 
3? {A;}. Furthermore, it is worth emphasizing that the 
distribution function HuMhi^) has three distinct power- 
law regions where the transition is given by the product 

riTa-. 



Q-1 

-1 



for flTa- < 1 

for ilT(j « 1 
for K r2r^. 



(5.15) 



We especially note that the high-frequency asymptote has 
fractal (self-similar) properties. Such fall-off also arises 
for Levy a-stable distributions. This might reveal infor- 
mation on the underlying physics. 

Keeping in mind that and To- may be regarded as 
break-frequencies, we note that fractional Zener time- 
parameters estimation don't really represent single relax- 
ation times (frequencies) as for discrete Debye models, 
but rather break-times (frequencies) around which the 
model characteristics change. 

An inversion recipe to find the analogy of Hyji^) given 
some attenuation law ak{u!) was presented in [88|, by ap- 
plication of an approach tightly related to . The small- 
attenuation assumption can at least for low frequencies 
often be reasonable for compressional wave propagation 
in biological tissue. By contrast, for shear-wave propa- 
gation in tissue, the attenuation is generally much more 
pronounced @. 



3. The a ^ /3 case 

For the more general situation when a ^ /3, inverse 
Fourier transform on both sides of (j5.12p instead gives 



H{t)Cn{K,{n)}{t) 











+ (io;)" J 
f (itj)"-( 






(iw)" 



it). 



(5.16) 



Subsequently exercising the inverse Fourier transforms 
gives 

H{t)Ko (-(t/r,)") 

- (rf /r;)t"-^+ii?„,„_^+i (-(t/T.)")l . (5.17) 



By recognizing in the equation above the Laplace trans- 
form relation (j7.26p of the Appendix, the distribution 
which we choose to denote k'i^ml(^) identified: 

= «o/a,i (f^,T-") - «:o(rf /t,")/,^„_^+i (17, r"") 

~ irn ' (T<,rj)2" -f- 2{T„n)°' cos(a7r) + l' 
• (r^rj)" sin(a7r) - (t,17)'^ sin(/37r) 



— '^'i/Ml(^)- 



(5.18) 

The fractional Zener wave equation (|4.7p may at a first 
glance hence be contained within the NSW framework 
of multiple relaxation Q , when assuming a continuum of 
relaxation mechanisms with the compressibility contribu- 
tion as given by the distribution (|5.18p . We 
note that as a consequence of the 6 < 1 criterion for the 
identity ((7?^ to be valid, the step from ((5Tf| to ([5TS)) 
is only correct for a < /?. 

On the other hand, for a continuous distribution of 
relaxation process contributions (|5.9p to be physically 
meaningful, the distribution Ki,{fl) must be non- negative 
for all 57. A closer investigation of the distribution 
'^Lml(^) (|5.18p above reveals that no matter how the 
non- negative parameters are set, it cannot be positive for 
all r2. This is in accordance with the a = (3 thermody- 
namical restriction discussed in Section PlI Bl 

This calls for a modification of the a (3 version of 
the fractional Zener wave equation (14.71) in order to make 
it thermodynamically admissible. Based on arguments 
presented in [68l |. we suggest to start out from a modi- 
fied, five-parameter form of the constitutive relation p.2p 
which is equivalent to an ansatz analyzed in ^90j and re- 
viewed in [48[: 



ait) + T, 



— En 



e(t) 



(5.19) 



where we have a < (3. A related model has been ap- 
plied to cell rheology [9l|. From the relation (|5.19p it is 
straightforward to construct a time-fractional wave equa- 
tion using the methodology of Section |IV] hence leading 
to 



c5 9t2 



at" 



^ dtp 



= 0. 
(5.20) 



We encourage researchers to execute the inverse trans- 
forms which yield the K^{il) NSW framework relaxation 
distribution corresponding to the above wave equation. 



4. Relation to the rheological relaxation time spectrum 

In viscoelasticity, a relaxation time spectrunij_7?(T), 
related to Kiy(51) is commonly studied (see e.g. [52j and 
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references therein). It is related to the complex modulus 
through 

/oo 
//(T)e-*/^dlnT. (5.21) 
-oo 

It may be shown that for the fractional Zener model when 
setting = r^^, the r-dependency of H{t) differs by a 
factor T from n^Mhi^) of (jS.lSp . Figs. 5 and 6 of [52] 
illustrate that a = f3 gives symmetric H{t), while a ^ p 
breaks the symmetry, most significantly far away from 
the peak region. 

B. Attenuation and phase velocity 

The decomposition of the frequency-dependent 
wavenumber into its real and imaginary parts, gives 
the phase velocity Cp{uj) = w/iR {^(lli)} and attenuation 
ak{uj) = — 5{fc(aj)}. In general, the attenuation and 
the phase velocity are thus given by insertion of the 
generalized compressibility into the dispersion relation 
as 

ak{i^) = -3{fc} = -w03^y3|^K(a;)| and 

CpH=c./5R{fc} = p-^/V5R{v^}- ^^'^^^ 

For the fractional Zener wave equation (j4.7p with 
a = /3, the attenuation expression (|5.22p can be com- 
bined with the fractional Zener compressibility hz{lo) of 
(j4.6p to get the attenuation. This results in three distinct 
frequency regimes of attenuation power-laws determined 
by the products To-il and r^ri Q: 

{^i+a low-frequencies, 
^i-o!/2 intermediate frequencies, (5.23) 
^i-Q high-frequencies. 

Below, the fractional Zener phase velocities and attenu- 
ations are further investigated numerically. The results 
from such calculations are compared to what is found 
by insertion of the distribution Ki/ml(^^) of (j5.14p into 
the NSW generalized compressibility integral formula 
(j5.9p . This generalized compressibility is finally applied 
to (j5.22p . from which afc(a;) and Cp{Lu) are found. 

We use the latter calculation method to explore the 
effect of letting the continuum of relaxation mechanisms 
populate only a bounded frequency interval, rather than 
the entire S [0,c»] region. 

Figure [5] compares attenuation curves, the frequency- 
dependent phase velocity, and the distribution KvM'Li^)- 
Note that for many applications, the ratio t„ /t^ is only 
slightly larger than one. This implies that the interme- 
diate regime becomes negligible. For attenuation in sea- 
water [92I and air [o^, one usually considers only three 
discrete relaxation contributions each with a = 1, which 
results in the familiar ak oc for low frequencies and 
constant attenuation for high frequencies. 




10"^ 10° 10^ ^(f 

COT 

a 




10"^ 10° 10^ 10° 

COT 

CT 




10"^ 10° 10^ 10° 

flT 



FIG. 2: Properties of the fractional Zener model exemplified 
for To- = IOOOtc and a = P = 0.5. Top pane: frequency- 
dependent attenuation as a function of normalized wave fre- 
quency uj ■ Ta, where the three power-law regimes are distin- 
guishable. Middle pane: Frequency-dependent phase velocity 
as a function of normalized wave frequency to ■ Tct- Bottom 
pane: the corresponding normalized effective compressibili- 
ties Ki,ml{^) of the continuum of relaxation processes as a 
function of normalized relaxation frequency £7 • To-. 



In the figures we observe that the parameter may be 
regarded as related to break frequencies, where the distri- 
bution of relaxation frequencies crosses over between dif- 
ferent power-law regimes of k^{^}). Notably the low- and 
high-frequency asymptotes of k^(J1), which both have 
fractal properties, are also well visible. The break fre- 
quencies also correspond to where the attenuation crosses 
over between different power-law regimes. 



C. Causality and finite phase speed 

In particular, we observe that according to the NSW 
paper Q any physical mechanism that fits into the mul- 
tiple relaxation framework corresponds to finite phase 
speed, non-negative attenuation, and causal response at 
all wave frequencies. This parallels the conversion of 
stress-strain models into the Maxwell- Wiechert represen- 
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tation. The NSW model causality is also verified because 
it complies with the Kramers-Kronig causality relations. 
See [9J] for an acoustics-oriented treatment of these rela- 
tions. 

In Section IV A 1[ we re-wrote the continuous relax- 
ation process distribution (j5.9p , which is a generalization 
of the NSW discrete sum (|5.8|) . into expression (|5.1ip : 

kn(w) = Ko - iujTt^H{t)Cn{K,^{^)}{t)^iLo). Studying 

this expression sparks the conclusion that any distribu- 
tion of relaxation contributions K^{n) for which the suc- 
cessive Laplace and Fourier integrals of (|5.11l) exist, the 
corresponding wave equation gives finite phase speed, 
non-negative attenuation, and causal solutions for all 
wave frequencies. Moreover, it is worth emphasizing that 
some models are causal but can imply unbounded phase 
speeds. For example the fractional Kelvin- Voigt wave 
equation the unbounded phase speed at infinite frequen- 
cies 0, m, HI] . The causality properties of several acous- 
tical attenuation models are investigated in [9]. Regard- 
ing restrictions on the attenuation power-law exponents, 
|96j argues that causality is maintained only if the at- 
tenuation has a slower than linear rise with frequency in 
the high-frequency limit. This requirement is met both 
for the fractional Zener wave attenuation and the NSW 
attenuation. See also [97] for a related study. 



erences therein. Maybe the search for first principles ex- 
planations for the fractional behavior of complex elastic 
media can be inspired by such findings. 

Furthermore, we note that relaxation processes in nu- 
clear magnetic resonance (NMR) as often described by 
the so called Bloch equations actually also can give 
non-Debye appearance, see e.g. jl02l ] and the references 
therein. In addition, developments related to dispersion 
and attenuation of elastic waves have many traits in com- 
mon with the mathematical descriptions of luminescence 
decay, see e.g. [103]. 

We aim to encourage the acoustical community to more 
frequently adopt fractional calculus descriptions for wave 
modeling in complex media. Hopefully we can also stim- 
ulate both mathematics-oriented and other researchers 
to spark further progress within fractional modeling of 
elastic waves by contributing to advance in e.g. model 
enhancements, existence and Green's function consider- 
ations, as well as in analytical and numerical investi- 
gations. We call for further exploration of connections 
between fractional dynamics, the surprisingly prevalent 
power-law patterns of nature, and the micromechanical 
structure of complex materials. 



VII. Appendices 



VI. Discussion and concluding remarks 

Within the framework of elastic waves, the current 
work surveys connections between concepts of power-law 
and multiple relaxation attenuation, causality, and frac- 
tional derivative differential equations. 

The fractional Zener elastic wave equation model fits 
attenuation measurements well and is characterized by 
a small number of parameters. The NSW model [3] is 
widely acknowledged in acoustical modeling and does not 
comprise fractional derivatives. Here we have analyzed 
how the fractional Zener (13. 2|) and NSW (15.91) wave equa- 
tion models can be unified under the assumption of a 
continuous distribution of relaxation mechanisms (|5.14p 
which has fractal properties. 

Because NSW-compatible wave equations are causal as 
well as consistent with non-negative attenuation and fi- 
nite phase speed, we prefer to consider any such model as 
eligible form an physically intuitive point of view. Never- 
theless, it is still unclear to the acoustics community what 
are the underlying physical mechanisms which interplay 
in complex media to result in power-law attenuation of 
elastic waves. 

The characteristics of viscoelasticity and acoustics of 
complex media, share many similarities with what is en- 
countered for dielectrical properties of complex media. 
We here point out that there are several theories within 
this field on how to explain e.g. Cole-Cole behavior by 
medium disordering, scaling, and geometry as well as 
more probabilistic approaches. See 98l - ll01 and the ref- 



A. Appendix: Mittag-Leffler function properties 

1. Definition and Fourier Transform Relation 

The on e-pa rameter Mittag-Leffler function was intro- 
duced in |104| . A two-parameter analogy was presented 
in jl05| . and may be written as 



E,,_h t] 



V{an -f h) 

n— 



(7.24) 



where F is the Euler Gamma function and the parameters 
are common ly re stricted to {a, 6} G C, 3?{a, 6} > 0, and 
i S C. See [lOGj ] for a comprehensive review of Mittag- 
Leffler function properties. 

A useful Fourie r tra nsform pair involving the Mittag- 
Leffler function is 

T{E[t) t'-'EaMi - Ar)} (i.) =^^^^. (7.25) 



2. Laplace Transform Integral Representation 

The function t^-'^Ea,b{-At'') ma y for < a < 6 < 1 
be written on an integral form jl08l ll09j |: 



6-1 



/>oo 

Ea,b{~Ae)^ e-"*/,.b(n,A)dl], (7.26) 
Jo 
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where 



Acknowledgements 



A sm[{b - a)7r] + sin(67r) 
TT + 2Art°- cos(a7r) + 



(7.27) 



Careful reading of Appendix E in [43| reveals that the 
above above functions fa,bi^, A) may be denoted spectral 
functions, which are non-negative. 
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